arXiv:1508.02376v2 [cond-mat.soft] 12 Aug 2015 


Brownian microhydrodynamics of active filaments 

Abhrajit Laskai^ and R. Adhikarjf] 

The Institute of Mathematical Sciences, CIT Campus, Chennai 600113, India 

Slender bodies capable of spontaneous motion in the absence of external actuation in an other¬ 
wise quiescent fluid are common in biological, physical and technological contexts. The interplay 
between the spontaneous fluid flow, Brownian motion, and the elasticity of the body presents a 
challenging fluid-structure interaction problem. Here, we model this problem by approximating the 
slender body as an elastic filament that can impose non-equilibrium velocities or stresses at the 
fluid-structure interface. We derive equations of motion for such an active filament by enforcing 
momentum conservation in the fluid-structure interaction and assuming slow viscous flow in the 
fluid. The fluid-structure interaction is obtained, to any desired degree of accuracy, through the so¬ 
lution of an integral equation. A simplified form of the equations of motion, that allows for efficient 
numerical solutions, is obtained by applying the Kirkwood-Riseman superposition approximation 
to the integral equation. We use this form of the equation of motion to study dynamical steady 
states in free and hinged minimally active filaments. Our model provides the foundation to study 
collective phenomena in momentum-conserving, Brownian, active filament suspensions. 


I. INTRODUCTION 

Slender bodies capable of spontaneous motion in vis¬ 
cous fluids are common in biological, chemical, physical 
and technological contexts. Examples from biology, in 
increasing degree of molecular complexity, include micro¬ 
tubules driven by molecular motors, axonemes, cilia, flag¬ 
ella mm, their synchronization and metachronal wave 
mm- In chemistry and physics, self-assembled bundles 
of microtubules driven by kinesin motors yields a model 
experimental system in which broken symmetry, collec¬ 
tive excitations, and topological defects can be studied 
out of equilibrium mm- In technology, much recent re¬ 
search has been directed towards the synthesis of slen¬ 
der bodies capable of spontaneous motion [T5HT71 . Such 
self-actuated slender bodies are expected to have many 
microfluidic and biomimetic applications [ 18] . 

Despite the great variety in both the structure of the 
body and its mechanism of self-actuation, the examples 
above have three features in common: the spontaneous 
motion of the slender body produces flow in the ambient 
fluid, the body is of a size sufficiently small to make Brow¬ 
nian fluctuations important, and the body resists defor¬ 
mation produced by the spontaneous flow and Brownian 
fluctuations. Any universal emergent behaviour in active 
slender bodies must appear from the interplay between 
fluid flow, Brownian fluctuations and the elasticity. Such 
systems present a new class of fluid-structure interaction 
problems. 

In this paper, we construct a theory of active slen¬ 
der bodies, by modeling them as filaments that enforce 
slip velocities or non-equilibrium stresses at the fluid- 
structure boundary. A multitude of microscopic mecha¬ 
nisms can produce such velocities or stresses. Our theory 
isolates the specific microscopic details of self-actuation 
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in the boundary conditions, from which universal, macro¬ 
scopic fluid flow can be generated. Such flow results from 
the exchange of momentum between the body and the 
fluid, and since no external forces act on the body or the 
fluid, the sum of their linear momenta is conserved. In 
the absence of external torques, the sum of their angular 
momenta is similarly conserved. These two constraints 
are explicitly taken into account when computing the 
fluid flow within our theory. 

The flow is computed through a discretization which 
replaces the continuous filament by a chain of spheres 
connected by non-linear springs. The spheres produce 
spontaneous hydrodynamic flow while the springs penal¬ 
ize changes in filament length and filament curvature. 
The antecedent of such a bead-spring discretization of a 
continuous filament traces back to Kramers, who used it 
to model the dynamics of a polymer. The crucial dif¬ 
ference between the model of Kramers and our adapta¬ 
tion of it is that the spheres in our theory produce spon¬ 
taneous flows. Each sphere must, therefore, be active. 
Such chains of active spheres have been used previously 
to model active filaments. In the earliest such model 
cm Ha, the spheres produce dipolar stresslet flows but 
are individually non-motile and are assumed so large that 
Brownian effects are negligible. In a subsequent contri¬ 
bution m, the spheres are taken to be motile, Brownian 
effects are included in two-dimensions, but contributions 
from non-local hydrodynamic flow are neglected. In a re¬ 
lated model, passive spheres are driven by tangential ac¬ 
tive stresses, hydrodynamically correlated Brownian mo¬ 
tion is included in three-dimensions, but active flow is 
neglected J22J [23] . Our theory presented here contains 
all previous models as special cases. 

In our recent work |24| , the problem of computing the 
fluid flow of N active bodies has been solved by trans¬ 
forming the local conservation law for momentum, which 
under the conditions relevant to the microhydrodynamic 
regime is the Stokes equation, into an integral equation 
over the sphere boundaries. The solution of this bound¬ 
ary integral equation gives the rigid body motion of the 
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active spheres, constrained by the conservation of linear 
and angular momentum, as a linear function of the forces, 
torques and the active boundary condition: 
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In the above, V„ and fl„ are the velocity and angular 
velocity of the n-th particle, V“ and are the self¬ 
propulsion and self-rotation contributions to the rigid 
body motion, the /i are the usual mobility matrices that 
relate external forces and torques to the rigid body mo¬ 
tion and the 7r are propulsion matrices |24| that relate 
Vm +1 \ the Z-th mode of the non-equilibrium slip veloc¬ 
ity on the surface of the m-th particle, to the rigid body 
motion (see below). These relations clearly show that 
rigid body motion of active particles is possible in the 
absence of external forces and torques, F( n = 0, T( n = 0, 
and even in the absence of self-propulsion, V“ = 0 and 
self-rotation = 0. Propulsion matrices, and not mo¬ 
bility matrices, are the key quantities that describe the 
correlated motion of active particles in a viscous fluid, 
constrained by the conservation of momentum and an¬ 
gular momentum [24]. While mobility and propulsion 
matrices have to be computed numerically for particles 
of general shape, analytical expressions can be derived 
when the particles are spheres [241 . 

Here, we use mobility and propulsion matrices for 
spheres of radius a, computed in the superposition ap¬ 
proximation first introduced by Kirkwood and Riseman 
[25], again, in the context of the dynamics of a poly¬ 
mer. In this approximation, the mobility matrices re¬ 
duce to the well-known Rotne-Prager-Yamakawa tensors 
while the propulsion matrices are obtained analytically 
[24] as gradients of the fundamental solution of the Stokes 
equation for an unbounded fluid. For a passive polymer, 
Yoshizaki and Yamakawa [2(3] verified that the superposi¬ 
tion approximation is correct to 0((a/b) 3 ), where b is the 
mean separation between spheres, as ./V —> oo . Since the 
propulsion matrices decay more rapidly with separation 
than mobility matrices, the superposition approximation 
for active filaments is also accurate to 0((a/b) 3 ) fM\ . 
The resulting equations of motion are used to study the 
dynamics of active filaments where the spheres produce 
dipolar flows and, so, are not individually motile. This 
part of our work complements studies of active filaments 
consisting of individually motile particles. We empha¬ 
sis that our general theory includes both the motile and 
non-motile cases studied previously. 


In the next section we present a generalization of the 
theory of Brownian motion of hydrodynamically inter¬ 
acting spherical particles to include surface activity. In 
section III, we construct equations of motion for active fil¬ 
aments using the results of this general theory. In section 
IV we introduce a minimal model for an active filament 
by discarding all but leading terms in the filament equa¬ 
tions of motion derived previously. Here we also study 
dynamics of such active filaments when both ends of the 
filament are free and when one end is tethered and the 
other end is free. In the first case, for sufficient strength 
of activity the filament is unstable to transverse pertur¬ 
bations, which results in a sequence of translational, ro¬ 
tational and oscillatory dynamical steady states. In the 
second case, we find a sequence of rotational and oscilla¬ 
tory dynamical steady states, both of which are limit cy¬ 
cles in the phase space of the dynamical system described 
by the equations of motion. In section V, linear stability 
analysis shows that the transition to a dynamic state hap¬ 
pens due a simple instability in both cases. This analysis 
shows that non-local hydrodynamic interactions are es¬ 
sential for the dynamic instability to occur. We conclude 
with a discussion of the application of our equations of 
motion to study collective phenomena in suspensions of 
active filaments [27; ;28| and other soft dissipative struc¬ 
tures. 


II. BROWNIAN MICROHYDRODYNAMICS OF 
ACTIVE SPHERES 

We consider the motion of N spherical active parti¬ 
cles of radius a in an incompressible viscous fluid. The 
center of the n-th sphere is at R n and its orientation 
is specified by the unit vector p„. The fluid can ex¬ 
change both momentum and angular momentum with 
the particles, of amounts determined by integrals of the 
momentum flux over the particle boundaries. In addi¬ 
tion to this contact contribution, the particles may be 
acted on by body forces and torques. At low Reynolds 
numbers inertia, of both the particles and the fluid, can 
be neglected, which results in an instantaneous balance 
of surface and body contributions to forces and torques. 
Newton’s laws, therefore, reduce to a pair of constraints 
on the fluid stress at the surface of every particle, 

MV„ = j> n • (TdS n + F e n = 0, (2a) 

m n = j> r x n • crdS n + T® = 0. (2b) 

The fluid stress er = — pI+77(Vv + Vv T ) has both hydro¬ 
static and viscous contributions and is determined from 
conditions of incompressibility and local momentum con¬ 
servation 


V • v = 0; V • a = 0. (3) 

where v is the fluid velocity, p is pressure and p is the 
viscosity. The solution of this Stokes system provides the 
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stress, from which the contact contribution of the force 
and torque on the every particle can be determined. In 
the absence of inertia and body forces, the conservation 
of particle momentum requires that the net contact force 
and the net contact torque on every particle be zero. 
The solution of the Stokes equation is determined by the 
boundary conditions on the particle surfaces and at in¬ 
finity. While activity can be expressed through a variety 
of boundary conditions on both the fluid velocity and on 
the fluid stress , we assume here an active slip at the sur¬ 
face of the particle [2H • This encompasses a wide variety 
of active phenomena, including electrophoresis, diffusio- 
phoresis [213, self-phoresis due to chemical catalysis m , 
and even swimming of microorganisms mm- We chose 
the fluid to be at rest at infinity. The boundary condi¬ 
tions, therefore, are 

v(R„ + p n ) = V n + r l n x p n + v a (p„) (4a) 
|v| —► 0, |p| 0, |r| —► oo. (4b) 

The task of obtaining the solution of the Stokes equa¬ 
tion is substantially simplified by recognizing that the 
three-dimensional partial differential equation can be re¬ 
duced, instead, to a two-dimensional integral equation 
over the particle boundaries. The starting point of this 
development is the integral representation of Stokes flow 

GUOS], 
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which provides the solution of Stokes equation in terms of 
the velocity and the traction, f = n • er, over the particle 
boundaries. The kernels in the integral representation 
are the Green’s function G, the pressure vector p and 
the stress tensor, K. 

Pi(r,r') = -V,V 2 p = ^ (6a) 

P 6 

Gij{ r,r') : (V% - VjVj) P = y + (6b) 

I^ijk (pi* ) — SikPj -f- V k G ij T Y iGjk ■ (6c) 

In the absence of boundaries, these kernels are transla- 
tionally invariant and, so, are functions of the separation 
p = r — r 7 . The boundary integrals can be evaluated an¬ 
alytically by expanding the boundary fields in complete 
orthogonal bases [21, 155], which are most conveniently 
chosen to be the tensorial spherical harmonics, 

Y«, 2 ... Ql (p) = (-1 ) /+ y +1 V Q1 ... V ai (£) . (7) 

In this basis, the surface velocity and traction on each 


particle is expanded as [551 55] 
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where the orthogonality of the basis functions yields the 
expansion coefficients as 

F - +1) = l\(2l — l)!! / f(Rm + Pm)YW(p m )dS m , (9a) 

V ™ +1) = Irm© / V(Rm + Pm)Y W (p m )dS m . (9b) 

These expansion coefficients are tensors of rank l + 1, 
symmetric and irreducible in their last l indices |24L ; 55| . 
The forces, torques, velocities and angular velocities are 
obtained from 

F£> = 2e • aFW = T e m 

= v m - <v“ >; 2e • v£> = -aVt m + ±(p m x <> 

where the angle brackets indicate integration of the en¬ 
closed term over the surface of the sphere and dividing 
the result by the surface area. Inserting the expansions in 
the boundary integral representation leads to a succinct 
expression for the fluid flow in terms of the expansion 
coefficients [21], 
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The boundary integrals G© and KW can be evaluated 
explicitly in terms of the Green’s function G and its 
derivatives. In this expression, the velocity coefficients 
can be computed from the boundary condition, but the 
traction coefficients remain unknown. To determine the 
traction coefficients, the boundary condition is first en¬ 
forced on the boundary of n-th particle, the resulting 
equation is weighted by the l-th tensorial harmonic and 
finally integrated over the n-th boundary. This Galerkin 
procedure yields an infinite-dimensional linear system of 
equations for the unknown traction coefficients [24], 
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where the matrix elements G[© and k(©,, 1 can, again, 
be evaluated analytically in terms of the Green’s function 
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G and its derivatives |M]. The formal solution of this lin¬ 
ear system yields Eq. ^ which expresses the rigid body 
motion in terms of the body forces and body torques and 
the coefficients of the active slip velocity. These equations 
reflect the linearity of both the Stokes equations and the 
boundary conditions. The hydrodynamically correlated 
motion is a sum of passive terms, proportional to the 
external body forces and torques, and active terms, pro¬ 
portional the modes of the active velocity. The last terms 
in each of the equations are the self-propulsion and self¬ 
rotation velocities of an isolated active sphere, which are 
given in terms of the active slip velocity as 

47ra 2 V“ = - J v a ( Pn )dS n (12a) 

47ra 2 0“ = -7^2 J Pn x v a (p n )dS n (12b) 

The above two equations have been known in the con¬ 
text of phoresis m and swimming by surface distortions 
[25]. They were later derived by the use of the Lorentz 
reciprocal relation [39]. The work of [35] derives these 
from the boundary integral representation of Stokes flow. 
The propulsion matrices emerge naturally from the solu¬ 
tion of the boundary integral equation as coefficients that 
determine the active contribution to the hydrodynamic 
interaction. Thus, Eq. [T] expresses hydrodynamic inter¬ 
actions between active particles acted upon by external 
forces and torques [24] . 

We note that the boundary integral method yields the 
hydrodynamic interaction between particles without the 
need to resolve bulk fluid degrees of freedom. This makes 


Eq. 1 especially useful for computing the hydrodynamic 
interaction of active particles in three dimensions, as the 
computational cost of resolving fluid degrees of freedom 
is completely eliminated [2T> IT) . 

The addition of Brownian fluctuations to Eq. |T] is ac¬ 
complished by appeal to linearity, the balance of dissipa¬ 
tion and fluctuation and Onsager symmetry of the mo¬ 
bilities. The generalization of the Einstein-Smoluchoswki 
description of the diffusion of a passive colloidal particle 
to the hydrodynamically correlated diffusive motion of N 
colloidal particles was completed by several authors us¬ 
ing Liouville, kinetic theory, Fokker-Planck and Langevin 
approaches HU- The Langevin approach provides the 
most direct way of incorporating in Eq. [T] by promoting 
them to a set of stochastic differential equations with a 
state-dependent fluctuation. The fluctuations are cho¬ 
sen to compensate each source of dissipation such that 
distribution of positions and orientations is Gibbsian in 
the absence of activity. The fluctuations, then, are cor¬ 
related Wiener processes with variances proportional to 
the mobility matrices. The positivity of dissipation en¬ 
sures that mobility matrices are positive-definite and On¬ 
sager symmetry constrains them to be symmetric in both 
the particle and translation-rotation indices. These two 
properties ensure that a mobility matrix can be factorised 
into a lower triangular matrix and its transpose, any one 
of which is a “square-root” of the mobility matrix. The 
fluctuations can then be expressed as products of uncor¬ 
related Wiener processes £ T , £ R , r] T , t) R and the “square- 
root” Cholesky factors. By linearity, dissipative, Brow¬ 
nian and active motions are additive. Therefore, the 
generalization of Brownian hydrodynamics to N active 
particles, expressed as coupled Langevin equations, is 
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These equations are the main result of this section. In 
the absence of activity, they reduce to the equation of 
Brownian dynamics with hydrodynamic interactions m- 
When the forces derive from positional and angular po¬ 
tentials, the form chosen for the fluctuations ensures that 
the Gibbs distribution of the positions and orientations 
is the stationary distribution. When activity is included, 
the balance between fluctuation and dissipation is no 
longer maintained and stationary states are no longer 
described by the Gibbs distribution. As we show in the 
remainder of the paper, non-trivial stationary states are 
obtained when the spheres are chained together into fil¬ 


aments. 


III. BROWNIAN MICROHYDRODYNAMICS 
OF ACTIVE FILAMENTS 

The equations of active Brownian hydrodynamics pre¬ 
sented in the previous section form the basis of our the¬ 
ory of active slender bodies. As mentioned before, we 
approximate the slender body as a filament, and then 
discretize the filament as a set of connected beads. The 
forces in the active Brownian hydrodynamic equations 
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Figure 1. Streamlines of irreducible modes of slow viscous flow around an active sphere. Displayed in (a) is • F, the flow 
around a sphere translating under the action of a force F; in (b) is J'VGOV' 2 ^ the minimal flow around a non-motile active 
sphere; and in (c) is V 2 G © the minimal flow around a motile active sphere. In the limit a — > 0, the flows in (a) and 

(b) reduce to a Stokeslet and stresslet, respectively, while (c) is a degenerate velocity quadrupole. The background colour is 
proportional to the logarithm of the velocity magnitude. 


then are derivatives of the non-linear spring potentials, 
F n = -V n U 

We assume the potential U to be a sum of connectivity, 
elastic, and self-avoidance potentials, 


N -1 N -1 

U = t/ C (R m , R m+ i) + P £ (Rm-l,Rra, Rm+l) 

m —1 m =2 

+ J2 US ( Rn,Rm) (14) 

m<n 

The connectivity potential U c (R, R ) = k{r - & 0 ) 2 /2, 
with elasticity parameter k, penalizes departures of the 
distance, r = |R — R |, of two consecutive particles from 
the equilibrium value of &o- The elastic potential for 
bending U E = R (1 — cos(f>) , with rigidity parameter /c, 
penalizes departures of the angle </> between consecutive 
bonds from their equilibrium value of zero. The rigidity 
parameter k connects to bending rigidity n as R = nbo. 
The self-avoidance potential U s restricts the overlap of 
particles and is taken here to be a Lennard-Jones po¬ 
tential that vanishes smoothly at a distance a^j = 2ser. 
The net force on the filament vanishes, as can be easily 
verified by summing the force on each particle. 

An approximate method for computing the hydrody¬ 
namic interactions, used widely in the dynamics of poly¬ 
mers, is to neglect their many-body character and, in¬ 
stead, assume them to be pair-wise additive. This su¬ 
perposition approximation, first introduced by Kirkwood 
and Riseman |25|, yields the well-known Rotne-Prager- 
Yamakawa form for the translational mobility. In that 
approximation, the mobility matrices of Eq. 1 are, 
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The diagonal parts of these matrices are one-body 
terms while the off-diagonal parts represent the hydrody¬ 
namic interactions. The diagonal parts are the familiar 
Stokes translational and rotational mobilities while the 
off-diagonal parts can be recognised as the Rotne-Prager- 
Yamakawa tensors and their generalizations to rotational 
motion. The Onsager symmetry of the mobility matrix is 
manifest in these expressions. In the same approximation 
the propulsion matrices are |24| 
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Figure 2. Streamlines of slow viscous flow around contractile (top row) and extensile (bottom row) minimally active filaments. 
Spontaneous flow in the linear conformation tends to compress contractile filaments and extend extensile filaments (first column). 
Both symmetric and antisymmetric transverse perturbations are suppressed in contractile filaments but enhanced in extensile 
filaments by the spontaneous flow (second and third columns). The latter leads to a linear instability in extensile filaments, 
when the elastic restoring force is no longer sufficient to counter the destabilizing tendency of the spontaneous flow. 


that accounts for the finite-size Faxen correction for the 
flow due to a sphere of radius a and c; is a numerical con¬ 
stant. The propulsion matrices, as defined here, represent 
hydrodynamic interactions and, therefore, are purely off- 
diagonal. The self-propulsion and self-rotation are wholly 
contained, respectively, in the diagonal terms V“ and 
Ct a . 

n 

For an active filament, the rotational motion of the 
spheres can be ignored, assuming torsional potentials 


that prevent such rotations. Therefore, the orientation 
of the spheres is no longer a dynamical degree of freedom 
to be determined from the angular velocity, but is pre¬ 
scribed. It is natural to fix the orientation of n-th sphere, 
p n = at„+/3n„+7b ra , in the local Frenet-Serret frame in 
terms of the of the tangent t„, normal n„ and binormal 
b„, and the direction cosines a , f3 and 7 . Combining all 
of the above, we obtain the following equation of motion 
for active filaments, 
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These stochastic differential equations describe the Brow- nian dynamics of an extensible, semi-flexible, self- 
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Figure 3. Non-equilibrium stationary states of a free, (a) - (f), and tethered, (g) - (i), minimally active filament. The figures 
show the filament conformation and the streamlines of slow viscous flow produced by the activity, where the background is 
coloured by the logarithm of the magnitude of the velocity. Free filament steady states are shown for three lengths with 
increasing values of activity _4 S . For short filaments (N = 48) increasing activity produces the three non-equilibrium steady 
states in (a) - (c), corresponding to the excitation of the first two elastic eigenmodes and their linear combinations. As the 
filament length is increased (N = 80 and N = 128) higher elastic eigenmodes appear with increasing amounts of activity in (d) 
- (f). Tethered filament steady states are shown for a fixed length (N = 64). Increasing activity produces a rotating steady 
state in (g), which bifurcates into an oscillating steady state in (h), with a return to a distinct rotating steady state in (i). 
(Movie) 


avoiding active filament including hydrodynamic interac¬ 
tions that arise from the exchange of momentum between 
the filament and its local conservation in the bulk fluid. 
These equations are the natural extension of Brownian 
hydrodynamics of passive filaments to the active case. In 
the next section, we study the simplest version of these 


equations, where only the most dominant contribution 
to active flow is retained and Brownian fluctuations are 
neglected. 


























































































IV. DYNAMICS OF MINIMALLY ACTIVE 
FILAMENTS 


In this section, we study in detail the simplest member 
of the family of models described by the general equa¬ 
tions of motion (Eq. 17) of the previous section. In this 


minimal model, the spheres are non-motile, V“ = 0, and 

all active velocity components other than the symmet- 
( 2 ) 

ric part of V„ are zero. Thus each individual sphere 
produces the flow show in Fig. lb. The velocities and 
tractions are, therefore, 


v (Rm T Pm) — R m T * pm 
47ra 2 f (R m + p m ) = -V m U + 3S, : 


(18a) 

(18b) 


where s m and S m are, respectively, the symmetric parts 
of the second-rank velocity and traction coefficients. The 
solution of the boundary integral equation, in the diag¬ 
onal approximation, relates the unknown traction coeffi¬ 
cient to the known value of the velocity coefficient [23] 05] 
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To complete the model, it is necessary to specify the ori¬ 
entation p n of the spheres, and hence the principal axis 
of the stresslet, in relation to the filament conformation. 
Motivated by the experimental observation that molecu¬ 
lar motors walking on microtubules generate tangential 
stresses [T], we parametrize s m uniaxially, with its prin¬ 
cipal axis always parallel to the local tangent t m of the 
filament, 


— So(t m t m 3 5). (^0) 

The coefficient So is positive for extensile stresses and 
negative for contractile stresses. Additionally, we assume 
that the activity is so large that the Brownian fluctu¬ 
ations make a negligible contribution to the dynamics. 
Active flow is balanced entirely by the filament elastic¬ 
ity. This leads to deterministic equations of motion for 
an active filament composed of non-motile beads, 


R„ = V T°T°G ■ V m U 
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These equations, without finite-sized corrections to the 
hydrodynamic flow, were first proposed in m and sub¬ 
sequently used in [2(J[ to study the dynamics of clamped 
active filaments. 

The flow produced by the filament is sum of contribu¬ 
tions from the potentials and the activity, 


N 


V r = 


v„i/. 


N 


7 4y.£I9°! 


Active 


( 22 ) 


The resultant flow is shown for three conformations, for 
both extensile and contractile filaments, in Fig 2. In the 
linear conformation, shown in the first column of Fig 2, 
the flow tends to compress contractile filaments and ex¬ 
tend extensile filaments. The stationary length of the fil¬ 
ament, then, is somewhat shorter in the contractile flow 
but somewhat longer in the extensile flow. In a symmetri¬ 
cally curved conformation, shown in the second column of 
Fig, 2, the spontaneous flow tends to suppress curvature 
in the contractile filaments but tends to enhance it for 
extensile filaments. The suppression and enhancement is 
seen for antisymmetrically curved conformations in the 
third row of Fig. 2. This shows that the interplay of flow 
and curvature is generally stabilizing for contractile fila¬ 
ments while it is destabilizing for extensile filaments. On 
dimensional grounds, a linear instability is expected when 
the filament length L > l a ~ k/sq. In the remainder of 
the paper, we shall focus only on extensile filaments and 
study the non-equilibrium stationary states that appear 
as a consequence of the linear instability. 

In d spatial dimensions, activity introduces a new rate 
T s = so/r]L d in addition to the rate of elastic relaxation 
= n/rjL d+1 of the filament bending modes. The ratio, 
A s = Lsq/k, of these two time scales provides a dimen¬ 
sionless measure of activity. The activity number A s also 
measures the departure from equilibrium and the amount 
of energy is injected into the fluid by the filament. We 
vary both the filament length and the activity number in 
studying the dynamics of the filament in d = 3 dimen¬ 
sions. 

We integrate the equations of filament equations of mo¬ 
tion numerically using a variable coefficient method. We 
chose the following parameters for the model : spring 
constant k = 1, bondlength b 0 = 4a, rigidity parameter 
R = 0.4, stresslet strength so = 0.0 — 0.5 and number 
of spheres N = 32 — 128. We obtain the mobility and 
propulsion matrices using the PyStokes library [4Q.I. We 
simulate the system for several hundred active relaxation 
time scales with an initial condition that is linear with 
small random, transverse perturbations. We study two 
cases, the first in which the filament is free at both ends 
and second in which it is free at one end and tethered at 
the other end. 

Our results are summarised in Fig. 3, which shows the 
non-equilibrium stationary states for both free filaments 
in panels (a) - (f) and for tethered filaments in panels (g) 
- (i). A linear instability appears at A s ~ 12 in free fil¬ 
aments which leads to spontaneous symmetric curvature 
and an emergent autonomous motility, shown in panel 
(a). The conformation corresponds, roughly, to the first 
elastic eigenmode of the passive filament. With increas¬ 
ing amounts of activity, higher elastic eigenmodes appear 
through a series of bifurcations, shown in panels (b) - (f). 
Whenever the conformation is asymmetrical about the 
center, the filament acquires a rotational component of 
motion. The higher elastic eigenmodes appear for smaller 
values of activity in longer filament, as is seen by com¬ 
paring panel (c) with panels (e) and (f). Although the 
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Figure 4. Comparison of filament kinematics and curvature dynamics. In (a) the ^-component of the center of mass velocity is 
compared with the ^-component of the scaled average curvature vector (see Eq. 23) for a free filament. In (b) the x-component 
of the mean angular velocity about the center of mas is compared with the scaled average first moment of the curvature vector 
(see Eq. 24). In (c) the same comparison is made for a tethered filament, with origin at the point of pivot. The comparison 
shows that the filament motion is strongly correlated with dynamics of activity induced curvature. 


system is three-dimensional, filament motion is planarly 
stable, in a plane that is determined by the initial condi¬ 
tion. This rich dynamics (movie I) and non-equilibrium 
steady sates we have found is remarkably similar to an 
experiment on an isolated axoneme [9]. 

Tethering the filament at one end, restricts transla¬ 
tion and thus the energy transduction from the activity 
is fed into rotational and oscillatory states (Movie II). 
Beyond the threshold activity of A s ~ 5, a tethered fila¬ 
ment rotates around the pivot, panel (g), in a plane that 
is chosen by the initial condition. Further energy injec¬ 
tion leads to flagella like beating in a plane, shown in 
panel (h). The highest value of activity studied is shown 
in panel (i), where a conformation similar to panel (d) 
appears, but now constrained by the pivot, is forced to 
rotate while maintaining conformation. With increase 
in filament length and activity, we expect higher elas¬ 
tic eigenmodes to appear, and either rotate or oscillate 
under the constraint of the tether. 

The dynamics of the center of mass follows from sum¬ 
ming the equation of motion over all spheres. The con¬ 
tribution from internal spring forces vanishes, and on ap¬ 
proximating the active flow by its contribution from the 
nearest neighbours, an approximate equation is obtained 
that relates the center of mass motion to the filament 


curvature and the activity, 

Vcm — (xn) = X (23) 

47T77O0 

In Fig. 4a, we compare the numerically computed values 
of the center of mass velocity with the curvature vector 
X defined above. There is a surprisingly good agreement 
between the two, indicating that principal effect of the 
non-local active flow can be expressed locally as a ten¬ 
dency to promote curvature. A similar relation holds for 
the angular velocity about the center of mass, 

flcM — - x (R — R cm)) = w (24) 

and the previous comparison is repeated for both free and 
tethered filaments in Fig 4b and Fig 4c. For tethered fila¬ 
ments the pivot point, and not the center of mass, is used 
to compute cross products, and the activity-dependent 
prefactor is two times smaller. This suggests that effec¬ 
tive, local equations of motion may be accurate for de¬ 
scribing some aspects of the dynamics of active filaments. 
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(a) Free filament with and without HI 



(b) Tethered Filament with and without HI 


Figure 5. Variation of the largest non-zero eigenvalues of the stability matrix, J, plotted against activity number, A s , for a 
free filament in (a) and a tethered filament in (b). In each case, the eigenvalues are computed including full hydrodynamic 
interactions (HI) (left panels) and neglecting all non-local hydrodynamic contributions (right panels). The eigenvalues remain 
negative, for all values of A s in a large range, when hydrodynamic interactions are neglected. Hydrodynamic interactions, 
therefore, are essential for the instability of the linear conformation and the bifurcation to dynamical steady states, signalled 
by the positive eigenvalues in each of the left panels. 


V. LINEAR STABILITY ANALYSIS dynamics, then, is 


To better understand the linear instability which is 
expected from the flows shown in Fig 2, we perform a 
stability analysis of the equations of motion, about the 
linear conformation. Taking the equations of motion to 
represent a dynamical system, R„ = /(R 1; R 2 , • • • , R n ), 


we compute the Jacobian J = 


-v n / 


r° 


ary state with linear conformation R°. 


at the station- 
The linearised 


5R n = J • JR„ (25) 

We numerically compute the eigenvalues of this stability 
matrix as a function of activity A s for both free and teth¬ 
ered filaments. To evaluate the importance of non-local 
hydrodynamic interactions, we also compare the eigen¬ 
values for the dynamics in which all non-local (that is 
m 7 ^ n) terms are deleted. The results are shown in 
Fig. 5(a) and Fig. 5(b) for free and tethered filaments, 
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respectively. 

We see that the largest eigenvalue becomes positive at 
A s ~ 12 for free filaments and A s ~ 5 for tethered fila¬ 
ments. The bifurcation is thus a simple instability and 
not of the Hopf type reported in the previous study of a 
clamped minimally active filament m- The first eigen¬ 
mode instability flows to the symmetrically curved con¬ 
formation shown in Fig 3 (a). Instabilities of the higher 
eigenmodes leads, the first of which is visible in Fig. 5a 
at A s ~ 40, produces the more complicated states shown 
in panels (b) - (g) of Fig 5. The presence of the tether 
increases the threshold value of the activity at which the 
instability occurs to A s ~ 5, but the sequence of insta¬ 
bilities remains identical. Remarkably, there is no insta¬ 
bility, in the same range of activity, when hydrodynamic 
interactions are ignored, as shown in the right panels of 
Fig 5. Thus, non-local active hydrodynamic flow is essen¬ 
tial to produce the instabilities and the non-equilibrium 
stationary states reported above. 

Subsequent bifurcations with increasing values of ac¬ 
tivity are expected to have a more complicated character, 
as the stationary states are typically limit cycles. The nu¬ 
merical study of limit cycle instabilities is considerably 
more involved than that of time-independent stationary 
states. We shall explore this aspect of the dynamics of 
active filaments in a future study. 


VI. DISCUSSION AND CONCLUSION 

Previous work on bead-spring models of active fila¬ 
ments have focused on three distinct mechanisms of ac¬ 
tivity. In the earliest work of Jayaraman et al m , ac¬ 
tivity arises from the hydrodynamic flow of the active 
spheres. In that work, the equations of motion for fila¬ 
ment dynamics in three dimensions contained contribu¬ 
tions from the leading order hydrodyamic flow due to 
stresslets and degenerate quadrupoles. A detailed study 
and all results were given for non-motile active spheres, 
thus ignoring the velocity quadrupoles. In subsequent 
work, Chellakot et al m studied a chain of motile active 
spheres, subject to Brownian motion, in two dimensions 
but ignored the all non-local hydrodynamic effects, both 
passive and active. In related work, Jiang and Hou J22] 
studied a chain of passive spheres, subject to forces of 
non-equilibrium origin, directed along the filament tan¬ 
gent. In their model, both passive hydrodynamic flow 
and hydrodynamically correlated Brownian motion is in¬ 
cluded in three dimensions, but active flow is absent. Re¬ 
markably, in spite of these differences between the mod¬ 
els, they yield a broadly similar phenomenology : linear 
instabilities, spontaneous motion, and oscillatory states. 

To understand why this might be, it is best to situate 
all the previous models within the equations of motion 
presented here. The model studied in detail by Jayara¬ 
man et al [19] is obtained when self-propulsion veloci¬ 
ties, V“, are set to zero, only the dipolar contribution 
to active flow is retained, and finite-sized corrections to 


hydrodynamic flow as well as Brownian motion are ne¬ 
glected. The model of Chellakot et al m is obtained 
when the self-propulsion velocity is directed along the 
axis p„, V“ = v s p n , and this axis is itself determined 
from the balance of a restoring and Brownian torques. 
This requires the angular velocity to be retained as a 
dynamical variable and all off-diagional contributions to 
mobility and propulsion matrices to be ignored. Finally, 
the model of Jiang and Hou m is obtained by ignoring 
all active components of flow, v a = 0, but representing 
the force on the spheres as F n = ~V„f7 + at n , as sum 
of contributions from the potentials and an unspecified 
non-equilibrium source. The common feature of all these 
models is that they produce motion in the direction of 
the curvature. This arises from the non-local hydrody¬ 
namic flow in the model of Jayaraman et al, and from 
the local contributions due to self-propulsion and non¬ 
equilibrium activity in the models of Chellakot et al and 
Jiang et al respectively. The present work shows that a 
phenomenology beyond curvature instabilities remains to 
be explored. In particular, torsional instabilities, possible 
with self-rotating active spheres that are unhindered by 
torsional potentials, are likely to yield further surprises 
in the dynamics of active filaments. 

The equations of motion presented provide the foun¬ 
dation for studying non-equilibrium statistical mechanics 
of active filaments. The coupled Langevin equation for 
the positions can be recast as Fokker-Planck equations, 
whose stationary solutions in the absence of activity are 
given by the Gibbs distribution. Activity, in the forms 
envisaged in this work, introduces an additional drift 
in the Fokker-Planck equation, destroying the balance 
between fluctuation and dissipation. This will lead to 
non-Gibbsian distributions in the stationary state, and, 
likely change well-known equilibrium properties like stat¬ 
ics of the coil-globule transition [43] and the distribution 
of loop closure times [30] ■ We urge that some of these 
problems be addressed both experimentally and through 
theory and simulations. 

As a final remark, we draw attention to the similarity 
between the instabilities reported here and the convective 
instability by active stress studied in the pioneering work 
of Finlayson and Scriven [41]. 
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Appendix A: Integrals and Matrix elements 

The expressions for the boundary integrals G'^ and and the matrix elements G nm } and K^nm are given as 


G^ +1 )( r , R m ) = J G(r,R m +p m )Y^(p m )dS m = a 1 2 * 4 AWjF z V«G(r, R m ) 
Kd +1 )( r , R m ) = n(2 ^ 1}!! J K(r, R m + p m ) • nY ^\p m )dS m = 


47ra a 1 

(l — 1)!(2Z + 1)!! 


(Ala) 

F l V«-VK(r,R m ) (Alb) 


G^+ 1 ’ z '+ 1 )(R„,R m ) = 


jf(j+i,i'+i)/p p» \ _ 

xv nra V xv n? ±x "mJ 


5 u 'liwa' / Y(0(/5) (^- pp ) y (/) ( p )^; 

J l+l ' 

—SwAttSA. 1 ' 1 '); 


fl“ l! '^>l i, V^G(R„,R„ l ); 


to = n; 

to 7 ^ n; 
to = n; 


4 7 ra p+r+ i ) 

(Z' - 1)!(2Z' + 1)!! 


V^'- 1 )K(R n , R m ); to ^ n; 


(A2a) 


(A2b) 


Appendix B: Constraining Force on the tethered bead 

Tethering the filament generates a conformation dependent constraint force at the tethered point. We take care of 
hydrodynamics to satisfy the boundary condition appropriately, such that net velocity of the respective end vanishes. 
Then we back-calculate the constrained force self-consistently on the course of simulation. 

1 1 N 7 c 3 N 

Rr = - (Vi U + F c ) - — Y, ^G(Ri,R m ) • V m U + — ^ VG(R l5 R m ) © s m 

m=2 m=2 

O JV TV 

F c = - Vi U - — ^°-A°G(Ri, R m ) • V m U + 7^ ^ J°J 1 VG(R 1 ,R,„) © s m . 

m—2 ra—2 
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